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Abstract: We report results to verify a theoretical framework to analyze the 
3D depth-wise structural organization of collagen fibers in articular cartilage 
using polarization-sensitive optical coherence tomography. Apparent 
birefringence data obtained from multi-angle measurements using a time 
domain polarization-sensitive optical coherence tomography system has 
been compared with simulated data based on the extended Jones matrix 
calculus. Experimental data has been shown to agree with the lamellar 
model previously proposed for the cartilage microstructure based on 
scanning electron microscopy data. This tool could have potential 
application in mapping the collagen structural orientation information of 
cartilage non-invasively during arthroscopy. 

© 2012 Optical Society of America 

OCIS codes: (110.4500) Optical coherence tomography; (170.3880) Medical and biological 
imaging; (1 10.5405) Polarimetric imaging. 
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1. Introduction 

Traditional techniques like arthroscopy, clinical magnetic resonance imaging (MRI), 
radiography and polarized light microscopy as diagnostic tools to study the onset and 
development of osteoarthritis have not been completely successful and each has its own 
drawbacks [1-3]. A non-invasive imaging tool to obtain information about the microstructure 
of the collagen fiber organization in articular cartilage has potential application in 
understanding and diagnosing different phases of osteoarthritis. Optical coherence 
tomography (OCT) is a promising candidate in this regard with its non-invasive ability and 
good resolution over penetration depths of l-2mm [4]. Specialized versions of OCT which 
make use of polarized light: polarization-sensitive optical coherence tomography (PS-OCT) 
offer more information on the structural orientation of the collagen fibers in articular cartilage 
[5,6]. Osteoarthritis is associated with the disruption of the delicate 3D structural organization 
of the Type II collagen fibers in the hyaline cartilage. Collagen fibers in articular cartilage 
tissues are organized in a complex 3D structural organization. A detailed study of the 
birefringence information obtained from light beam incident on the sample at multiple angles 
of incidence provides a more complete picture of the 3D orientation of the collagen fiber fast 
axis. This has been first shown by Ugryumova et al, for a single layered model of collagen 
fibers in tendon tissues [7]. A more rigorous theoretical basis to study the propagation of the 
polarized components of light in layered anisotropic media is based on the extended Jones 
matrix calculus (EJMC), developed by Yeh et al, for describing light propagation in liquid 
crystals. The EJMC describes light propagation through layered birefringent media with an 
arbitrary 3-D orientation of the optic axis in each layer by extending the conventional 2x2 
Jones calculus to the case of off axis light transmission, accounting also for the Fresnel 
transmission coefficients at the initial interface [8,9]. This provides a framework for studying 
layered polarized light propagation in articular cartilage, which possesses a complex 3D 
orientation of collagen fibers. The EJMC has been used by our group previously to 
theoretically model light transport in a scattering medium using a Monte Carlo method [10]. 
Also, a purely theoretical description pertaining to cartilage using the EJMC has been reported 
recently by Fanjul-Velez et al. [11]. They have also applied this approach to carry out a 
polarimetric study of the human cornea over its layered structure. 

In this paper, angle- and depth-resolved retardance profiles obtained from bovine articular 
cartilage using a time-domain PS-OCT system are interpreted using a layered model, which is 
based on scanning electron microscopy (SEM) data [12] and then propagated through the 
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EJMC. The 3-D geometry of the collagen fibers is then crudely estimated by parameterizing it 
and then manually adjusting these parameters to obtain the least-squares best fit of the data. 
To the best of our knowledge, this is for the first time a comparative study based on 
experimental and theoretical data of microstructure of articular cartilage has been reported in 
order to explain the PS -OCT data obtained. 

▲ 




Fig. 1 . (a) A schematic showing the EJMC approach to modeling the zonal layered structure of 
articular cartilage. Here (4 and 8 C give the azimuthal and the polar angle of the collagen fiber 
direction, i.e., the optic axis. The incident angle 6>, m - is varied with respect to the surface normal 
of the tissue over two orthogonal planes to obtain the multi-angle PS-OCT data. 1(b) A 
schematic of the imaging procedure of multi-angle PS-OCT study of articular cartilage over 
two orthogonal planes. 

2. Materials and Methods 

The layered morphology of the collagen fiber network in articular cartilage is as shown in Fig. 
1(a) and also shown in more detail in Fig. 2(b). The surface 'tangential' layer comprises 
typically 10% of the total thickness of the articular cartilage and consists of collagen fibers 
oriented parallel to the surface. This layer is followed by the 'transitional' zone, which as per 
the lamellar model of Jeffery et al. has collagen fibers organized into parallel sheets which 
arch downwards leading to the radial zone, in which the collagen fibers are oriented 
perpendicular to the surface [12]. The transitional zone comprises 40-60% of the total 
thickness of the cartilage with the rest of the thickness comprising of the radial zone leading to 
the subchondral bone [13]. We base our theoretical model of the articular cartilage on this 
leaf-like lamellar model of Jeffery et al., starting with a simple model with the assumption of 
constant azimuth angle orientation of the collagen fibers throughout the depth of the tissue 
[12]. Fresh tissue samples of bovine articular cartilage were extracted from the fetlock joint of 
the hindlimb of the animal obtained from the local abattoir. The samples were stored frozen at 
-20°C prior to imaging and then sectioned along the anterior side of the apex for imaging. 
Multi-angle (0°, +60° and -60° relative to the surface normal) PS-OCT measurements were 
carried out for two different orthogonal planes, x-zand y — z (Fig. 1(b)). The 
x, y, z Cartesian coordinate system is defined such that the z-axis represents the axial depth 
direction. The x— z and y- z planes corresponds to the coronal and sagittal planes, defined as 
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the plane in which the two connected bones are constrained to move and the orthogonal plane, 
respectively. Hence in general, these planes are oriented at an unknown angle with respect to 
the azimuthal orientation of the superficial collagen fibers. 
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Fig. 2. (a) A schematic of the bulk-optics TD-PS-OCT system used to carry out the study; 
QWP, quarter wave plate, BS, beam splitter, PBS, polarizing beam splitter, P, polarizer, (b) A 
schematic of the cartilage zonal layered structure and the layer thickness used for the EJMC 
study. Also shown are the orientations of the polar angle of the collagen fast axis varying from 
90° in the superficial zone to gradually becoming 0° in the radial zone. 

The theory of the time-domain PS-OCT (TD-PS-OCT) system used for our study is based, 
is similar to the first PS-OCT system demonstrated by Hee et al. [5], and other first generation 
TD-PS-OCT systems reported subsequently [14]. Our system is described in details elsewhere 
[15]. A schematic of the experimental set up is as shown in Fig. 2(a). The bulk optics TD-PS- 
OCT system uses a light source based on quantum dot superluminescent diode centered at 
1300nm with a bandwidth of 85nm [16]. The light beam is linearly polarized and then split 
into the reference arm and the sample arm. In the reference arm, a quarter wave plate is placed 
with fast axis oriented at 22.5°, which on forward and backward propagation along the 
reference path yields orthogonal components of equal amplitude. In the sample arm, the 
linearly polarized light is incident on a quarter wave plate oriented at 45° yielding light 
incident on the sample that is circularly polarized. The Fourier delay line scanning speed is 
100Hz which means it takes 4s to acquire an OCT image of 400 lateral scans (A-scans) over a 
width of 4mm. From the intensities detected at the two detectors placed after a polarizing 
beam splitter, I v (z) and I H (z) , the depth dependent retardance values are obtained as 



8{z) = arctan 



Iyiz) 

I B (z> 



(1) 



The nondepolarizing nature of the backscattering measurement involved in OCT allows 
Jones matrix based analysis of the measured signals [17]. Henceforth, we base our simulation 
on Jones matrix analysis, incorporating the extended Jones matrix calculus which takes into 
account off axis light propagation in birefringent multi-layered structures. We can model the 
light beam path of the circularly polarized light incident on the sample and then back onto the 
detectors through Jones calculus as 



A H (z) 
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where A H (z) and A v (z) are the electric field horizontal and vertical amplitudes of the light 
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beam on the detectors. QWP represents the quarter wave plate inserted in the sample. J k is 

the round-trip Jones matrix for the layered cartilage model, which we calculate using the 
EJMC. EJMC can be applied to the birefringent-layered cartilage model on the assumption of 
weak birefringence ( An = \n e -n o \ « n o ,n e ), which is satisfied for biological tissues such as 

cartilage. n e and n o are the refractive indices of extraordinary and ordinary modes of light 
propagation, respectively. Here, we present the analysis based on the electric field 
components of the two orthogonal components of the light beam detected. Although in 
principle, phase-resolved data analysis can yield complex depth-resolved amplitudes, in 
practice, the depth dependent retardance profiles are calculated using only the modulus value 
of the amplitudes, which then restricts the values of the retardances calculated to the range 
[0 tt/2\ for the experimental data. 

With the assumption of negligible multiple reflection occurring at the interface of the 
different layers in the cartilage model, the overall layered cartilage model is given by 



J^e^P'PT. (3) 



where 7] = 
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and T = 
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are the Fresnel reflection coefficients at the air tissue 

pj 

interface and P is the single-pass cumulative Jones matrix of the cartilage from the surface to 
the desired depth. For a structure comprising of m parallel layers, P can be represented as a 
product of Jones matrices of individual layers with each layer treated as a linear retarder with 
a fixed fast axis orientation: 
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coordinate system from that describing the normal modes of propagation (i.e., the o-wave and 
e-wave) to that describing the s-wave and p-wave. The full equations for T i ,T o ,R{y/ i ),k oz and 
k a can be found in [8] and [9]. In our model, the layered structure of cartilage is described by 
dividing the full cartilage thickness into m = 22 layers each of thickness d j (Fig. 2(b)). The 

superficial layer is assumed to be lOOum thick divided into 4 layers of uniform thickness. The 
transitional and radial layers are 500|im and 400um respectively with each further divided into 
layers of thickness 50|im. The polar angle of the superficial zone is 90° which gradually 
changes in the transitional zonal layers becoming 0° in the radial layer. The polar angle 
change across the depth of the tissue is based on the SEM studies [12]. In this study, the 
values of the polar angle input for the transitional and radial zone (for the odd numbered 
layers) are [85, 60, 45, 35, 28, 15, 0, 0, 0] degrees. All layers are assigned the same azimuthal 
fiber orientation along the depth of the tissue, an arrangement suggested by the lamellar 
arrangement of the fibers shown in the SEM images. A small azimuthal angle in the range of 
0-10° seems to provide a closer fit to the PS -OCT data than zero (azimuthal angle of 0° was 
used for this simulation); recall that the azimuth is unknown a priori as it is well known that 
the orientation of split-lines in articular cartilage is not necessarily parallel to the sagittal 
ridge. The refractive indices of the ordinary and extraordinary modes of light beam 
propagation are assumed to be constant over the depth of the cartilage. Whilst this assumption 
could be questioned, our aim in this study is a preliminary investigation to establish the 
applicability of this type of modeling to articular cartilage and many refinements of the model 
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are possible due to its inherent flexibility. The values of birefringence reported in the literature 
varies from (0.45±0.12)xlCT 3 for 4 month old porcine sample to (1.48±0.55)xl0~ 3 for a 21 

month old porcine sample [18] as well as a value of 3 xl0~ 3 reported by Jiao et ah, for bovine 
sample [19]. The theoretical model previously applied by Fanjul-Velez et al, assumes a high 
value of birefringence ( ~ 5xlCT 3 ) for articular cartilage, which is much higher than the value 
we find for bovine cartilage experimentally [11]. A value of An = 2.2 x 10~ 3 gave a good 
manual fit to the experimental data for the other parameters as stated. The polar angle of the 
fibers varies in the manner described above. All these parameters are propagated through Eqs. 
(2)-(4) to derive A H (z) and A v (z) . From these quantities the single-pass phase retardance in 
the noise-free case could be calculated using Eq. (5) [5]: 



S(z) = tan" 



(5) 



In practice, however, Eq. (5) cannot fully describe the experimental data because of the 
effects of measurement noise. A noise model based on the OCT signal attenuation in the 
single scattering approximation I(z) oc exp(-2/v ( z) is introduced. Here, /v ( is the attenuation 

coefficient. A noise bias term a is also introduced into the noise model similar to the 
background noise reported by Schoenenberger et al. [20]. The modified equation for the 
depth-dependent retardance is then given as 



S(z) = tan 



^[Mzf^i-l^ + a 2 ] 
^[A H (z) 2 cxp(-2 Ml z) + (J 2 ] 



(6) 



3. Results and Discussions 



For the above mentioned choice of parameters, angle-resolved depth-dependent retardance 
profiles for oblique incidences of ± 60° as well as for normal incidence in two orthogonal 
planes in the chosen coordinate system are simulated using the EJMC and closeness to the 
experimental data is visually assessed to obtain a good manual fit. This procedure of manual 
fitting was chosen initially to assess the likely success of an automated fitting algorithm based 
on nonlinear optimization. 
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Fig. 3. (a)-(e) Depth-dependent retardance profiles obtained from PS-OCT measurements 
(dots) compared with those obtained from EJMC simulation results (line + star) over multiple 
angles of illumination in sagittal and coronal planes of the bovine cartilage sample using 
manual fit. 

The simulated profiles are compared with the experimental data in Figs. 3(a)-(e) as 
obtained by fitting the parameters described and visually assessed for 5measurements of 
angle-resolved PS-OCT measurements together. Note that the measured profiles are obtained 
by averaging 50 neighboring A-scans to reduce the speckle noise. The range of the retardance 

71 



value is 



0 



as required by the inverse tangent function in Eq. (5) and the positivity 



ofA^(z) and A H (z). Wrapping of values that lie outside this range creates the obtained 

retardance profiles. Inclusion of the noise model as given in Eq. (6) makes the retardance 
profiles wrap before reaching the extremes of this range since neither the numerator nor the 
denominator in Eq. (6) can reach zero in the presence of noise bias. As plotted in Fig. 3(c), for 
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normal incidence as the incident beam direction is perpendicular to the collagen fiber optic 
axis direction in the superficial zone but parallel in the deeper radial zone, the cumulative 
retardance tends to display a strong band near the surface and remain constant in the deeper 
region, an observation which is supported by EJMC simulation. Figures 3(a)-(e) reveals how 
the banding pattern for retardance profiles over oblique angle of incidence along the two 
orthogonal planes of imaging obtained from PS-OCT experimental data tend to agree with the 
results obtained from EJMC simulation for our assumed cartilage model. A reasonable manual 
fit is obtained between the predicted retardance profiles and the experimental data. Note how 
the plane in which sheets are folded (i.e. the plane containing the surface normal of the sheets) 
can be located since incident illumination directions lying in this plane yield the maximum 
difference in banding periodicity (Figs. 3(a) and 3(b)). The result presented here has to be 
further supported based on histology studies to fully validate this technique as a tool to study 
the microstructure of articular cartilage. This work is ongoing in our lab. 

To extend the technique towards potential clinical applications an automated fitting 
procedure has also been developed using the simplex optimization algorithm 'fminsearch' 
from optimtool toolbox of Matlab™. It is based on Nelder-Mead direct search optimization 
method which works by minimizing an objective function [21]. The routine optimizes 8 free 
parameters: the thickness of superficial zone, t , the thickness of transitional zone, t tmm , the 
quadratic and linear coefficients 'a\ V of a quadratic function to describe the arching of the 
collagen fibers in the transitional zone, the true birefringence, n e - n o , the azimuthal 

orientation of the optic axis, <j> c , the tissue attenuation coefficient, /u t , and the noise bias term, 
a . The error metric (Eq. (7)) that is minimized is 

5 n 2 

X = 2j 2j [^'.measured i Z i) ~ ^/.modelled { Z i )] 0) 
j=l i=l 

where n denotes the number of axial measurement points, and j labels the incident k- 
vectors. The depth dependent change of 8 c is modeled as a quadratic profile (Eq. (8)) in the 
transitional zone as 

^an S (z) = a-z 2 +b-z + c (8) 

where z represents the depth into the transitional layer. The initial values of a, b, and c are 
obtained from the tabulated 6 C values used in the manual fit described earlier. Also, the value 



Table 1. List of all the parameters used in the 'fminsearch' optimizer for obtaining a good 
fit of the depth-dependent retardance profiles obtained from angle-resolved PS-OCT 
imaging and those obtained from EJMC simulation" 



Parameters 


Input condition 


Output obtained 


Thickness of superficial layer (in urn) 


50 


83 ± 27 


Thickness of transitional layer (in urn) 


600 


595 ± 360 


Azimuthal angle ((& in deg.) 


-25 


2.2 ± 8.2 


Coefficient a 


140 x 10" 6 


(115 ± 62) x 10" 6 


Coefficient b 


-21.77 x 10~ 2 


(-19.25 ± 10.46) x 10~ 2 


True birefringence An 


1.5 x 10 3 


(2.18 ± 0.067) x 10~ 3 


Attenuation coefficient, jj, 


4.0 x 10~ 3 


(2.5 ± 0.05) x 10~ 3 


Noise bias, a 


6.0 x 10~ 2 


(8.3 ± 4.1) x 10~ 2 



"This optimization would thus yield the depth-dependent 3D structural information data from the articular cartilage as 
shown in the Output column of the table. Also given are the values of the input condition given to the optimizer. 
Output 1 is obtained by averaging the output obtained by running the optimizer through 100 random generated data 
set derived from the averaged A-scan of PS-OCT using bootstrap Monte Carlo method. 
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Fig. 4. (a)-(e) Depth-dependent retardance profiles obtained from PS-OCT measurements 
(dots) compared with those obtained from EJMC simulation results (line + star) over multiple 
angles of illumination in sagittal and coronal planes of the bovine cartilage sample based on 
parameters obtained from 'fminsearch' optimizer and bootstrap Monte Carlo method. 



of 9 c gradually varies from 90° to radial along the transitional layer which restricts the 

coefficient 'c' to small bound [85 90] which typically gives R value at the start of the 

transitional layer. However, in this study, coefficient 'c ' in the quadratic profile of 0 is not a 

variable input into the optimizer, but is input as a constant of value 90 which describes 0 c at 

the onset of the transitional zone. The list of the parameters used in the optimizer and their 
initial and final values are tabulated in Table 1 and shown in Figs. 4(a)-(e). A noticeable 
improvement in the fit is produced by the optimizer although significant differences still 
persist and indicate that the model is still in need of further refinement. 
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The optimizer can also be used to estimate the precision on the fitted parameters, using the 
"Bootstrap Monte Carlo" procedure. We apply the optimizer to an ensemble of 100 random 
data sets, generated using the mean experimental data set and the estimated standard 
deviations at each axial point [22]. The statistics of the fitted parameters values are then 
measured directly. The mean profiles are averaged over 50 neighboring A-scans. To estimate 
the standard deviations we extract 10 A-scans randomly from the original set of 50 and 
generate one realization of the average. This is then repeated 10 times and the standard 
deviation at each axial pixel measured. The standard deviation for the 50 A-scan averages is 
then taken to be 1/V5 of this value. The resulting mean values along with the standard 
deviations of the 9 fitted parameters are tabulated in Table 1 . It is also noted that, within the 
set of constraints to the parameters given by physically possible values, a reasonable choice of 
any input conditions yield the output within the error range presented. 

In summary, we have shown the applicability of the EJMC approach to understanding the 
complex 3D structure possessed by articular cartilage by analyzing the PS-OCT data obtained 
for multiple angles of illumination. Although the PS-OCT system that we have used is not 
state of the art, the analysis shown in this paper is quite general. Further studies are underway 
to use this EJMC approach to quantitatively map the topographical variations in the 3D 
structural over the entire articulating surface, which ultimately could help in better 
understanding the biomechanics of the tissue and also detecting pathological changes. We 
have also applied this EJMC simulation to PS-OCT data obtained from a recently 
implemented continuous polarization modulation swept source PS-OCT system, the speed of 
which could potentially allow these measurements to be made in vivo. Results will be reported 
in due course. 
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